0. Setup

You’ll need to download Stan to your machine for these analyses. Note: this document will save brms model fits into your working directory.

First, let’s load all the packages we’ll need.

library(bayesplot)
library(brms)
library(cowplot)
library(ggplot2)
library(grid)
library(tidyverse)

Now, read in the data from Study 2.

(
  d <- 
    read.csv('data/exp2data.csv', header = TRUE) %>%
    as_tibble()
)

Each row is a different participant in Study 2 (there were 80). Columns 1-68 are variables that summarise the whole session. Columns 69-493 refer to events in each round of the Visible condition. Columns 494-918 refer to the Hidden condition. This experiment followed a within-subjects design (i.e. participants completed both the Hidden and Visible conditions). In the counterbalancing column, 0 means that the participant played the Visible game before the Hidden game, and 1 means that the participant played the Hidden game before the Visible game.

1. Binomial - Number of shocks

Trim the dataset and get it in long-format.

(
  d.trim <-
    d %>%
    select(ID, ends_with('.rounds_survived'), ends_with('.overall_num_shocks')) %>%
    gather(key, value, -ID) %>%
    extract(key, c("Condition", "variable"), 
            "(hidden|visible).(rounds_survived|overall_num_shocks)") %>%
    spread(variable, value) %>%
    mutate(Condition = ifelse(Condition == 'hidden', 1, 0))
)

As in Study 1, we first test whether the probability of a shock occurring differs between the two conditions. Although this experiment followed a within-subjects design, random effects are not necessary here because the occurrence of shocks are independent between the two conditions.

All analyses in this document will use the brm() function. Learn more about the brms package here. Also, see here for an explanation of why I use the 0 + intercept syntax throughout this document.

b2.1 <-
  brm(data = d.trim, family = binomial,
      overall_num_shocks | trials(rounds_survived) ~ 0 + intercept + Condition,
      prior = c(prior(normal(0, 1), class = b)),
      iter = 2000, warmup = 1000, chains = 4, cores = 4,
      seed = 2113)

save(b2.1, file = "models/exp2_brmsfit1.rda")

We set the seed to a random number, to make the results reproducible. Here are the priors that were set for this model.

prior_summary(b2.1)

Let’s look at the results.

print(b2.1)
##  Family: binomial 
##   Links: mu = logit 
## Formula: overall_num_shocks | trials(rounds_survived) ~ 0 + intercept + Condition 
##    Data: d.trim (Number of observations: 160) 
## Samples: 4 chains, each with iter = 2000; warmup = 1000; thin = 1;
##          total post-warmup samples = 4000
## 
## Population-Level Effects: 
##           Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## intercept    -1.36      0.06    -1.47    -1.25 1.00     1223     1649
## Condition    -0.03      0.09    -0.21     0.14 1.00     1244     1612
## 
## Samples were drawn using sampling(NUTS). For each parameter, Eff.Sample 
## is a crude measure of effective sample size, and Rhat is the potential 
## scale reduction factor on split chains (at convergence, Rhat = 1).

Plotting the parameters.

stanplot(b2.1)

The Condition parameter is -0.03, and its 95% credible intervals cross 0, implying that the probability of a shock is no different across the two conditions. On the probability scale:

post <- posterior_samples(b2.1)

visible_prob <- inv_logit_scaled(post$b_intercept)
hidden_prob  <- inv_logit_scaled(post$b_intercept + post$b_Condition)
difference   <- hidden_prob - visible_prob

quantile(difference,c(0.025,0.5,0.975)) %>% round(2)
##  2.5%   50% 97.5% 
## -0.03 -0.01  0.02

2. Gaussian - Total amount of cattle lost due to shocks

Trim the dataset again.

(
  d.trim <-
    d %>%
    select(hidden.total_cattle_lost, visible.total_cattle_lost) %>%
    gather(Condition, total_cattle_lost) %>%
    mutate(Condition = ifelse(Condition == 'hidden.total_cattle_lost', 1, 0))
)

We now fit a model to determine if the total amount of cattle lost due to shocks varies between conditions. Again, no random effects are needed because the outcome is independent across conditions (generated stochastically by the game), despite the within-subject design.

b2.2 <-
  brm(data = d.trim, family = gaussian,
      total_cattle_lost ~ 0 + intercept + Condition,
      prior = c(prior(normal(0, 100), class = b, coef = 'intercept'),
                prior(normal(0, 5), class = b, coef = 'Condition')),
      iter = 2000, warmup = 1000, chains = 4, cores = 4,
      seed = 2113)

save(b2.2, file = "models/exp2_brmsfit2.rda")

The priors we used.

prior_summary(b2.2)

The results.

print(b2.2)
##  Family: gaussian 
##   Links: mu = identity; sigma = identity 
## Formula: total_cattle_lost ~ 0 + intercept + Condition 
##    Data: d.trim (Number of observations: 160) 
## Samples: 4 chains, each with iter = 2000; warmup = 1000; thin = 1;
##          total post-warmup samples = 4000
## 
## Population-Level Effects: 
##           Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## intercept    51.27      2.11    47.26    55.53 1.00     2638     2655
## Condition    -6.70      2.70   -11.98    -1.30 1.00     2736     2727
## 
## Family Specific Parameters: 
##       Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## sigma    20.41      1.16    18.30    22.78 1.00     2617     2448
## 
## Samples were drawn using sampling(NUTS). For each parameter, Eff.Sample 
## is a crude measure of effective sample size, and Rhat is the potential 
## scale reduction factor on split chains (at convergence, Rhat = 1).

Plot the parameters.

stanplot(b2.2)

On average, 7 more cattle were lost in the Hidden condition than the Visible condition.

3. Binomial - Probability of requesting

Get a long-format data frame with binary request decisions over all rounds, for both conditions. If request == NA, player has died and been removed from the game, so we drop those rows.

(
  d.trim <-
    d %>%
    select(ID, Group, ends_with('.player.request')) %>%
    gather(key, request, -ID, -Group) %>%
    extract(key, c("Condition", "round_number"), 
            "SurvivalGame_(Hidden|Visible)[.](.|..)[.]player[.]request") %>%
    mutate(Condition    = ifelse(Condition == 'Hidden', 1, 0),
           round_number = as.integer(round_number)) %>%
    drop_na()
)

This leaves us with 2834 request decisions.

We now fit a varying intercept and slope model, with participants nested within groups. We allow the slopes for both round number and condition to vary, to fit with the experiment’s within-subjects design (participants completed multiple rounds, in both conditions).

b2.3 <-
  brm(data = d.trim, family = bernoulli,
      request ~ 0 + intercept + round_number + Condition
      + (0 + intercept + round_number + Condition | Group/ID),
      prior = c(prior(normal(0, 1), class = b),
                prior(student_t(3, 0, 10), class = sd),
                prior(lkj(1), class = cor)),
      iter = 2000, warmup = 1000, chains = 4, cores = 4,
      sample_prior = TRUE,
      control = list(adapt_delta = 0.99),
      seed = 2113)

save(b2.3, file = "models/exp2_brmsfit3.rda")

Here are the priors for the model we just fitted.

prior_summary(b2.3)

Now let’s see the results.

print(b2.3)
##  Family: bernoulli 
##   Links: mu = logit 
## Formula: request ~ 0 + intercept + round_number + Condition + (0 + intercept + round_number + Condition | Group/ID) 
##    Data: d.trim (Number of observations: 2834) 
## Samples: 4 chains, each with iter = 2000; warmup = 1000; thin = 1;
##          total post-warmup samples = 4000
## 
## Group-Level Effects: 
## ~Group (Number of levels: 40) 
##                             Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## sd(intercept)                   0.65      0.29     0.06     1.17 1.02      288      278
## sd(round_number)                0.05      0.02     0.01     0.08 1.01      387      502
## sd(Condition)                   1.65      0.31     1.07     2.30 1.00      939     1066
## cor(intercept,round_number)    -0.40      0.38    -0.89     0.56 1.01      441      918
## cor(intercept,Condition)       -0.25      0.33    -0.79     0.53 1.01      330      265
## cor(round_number,Condition)     0.34      0.30    -0.35     0.85 1.01      378      476
## 
## ~Group:ID (Number of levels: 80) 
##                             Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## sd(intercept)                   0.68      0.26     0.16     1.19 1.02      293      509
## sd(round_number)                0.04      0.02     0.00     0.08 1.02      251      779
## sd(Condition)                   1.01      0.32     0.43     1.69 1.01      432      686
## cor(intercept,round_number)    -0.41      0.44    -0.92     0.70 1.00      498     1404
## cor(intercept,Condition)       -0.32      0.34    -0.83     0.49 1.01      509      630
## cor(round_number,Condition)     0.11      0.44    -0.78     0.84 1.01      315      602
## 
## Population-Level Effects: 
##              Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## intercept       -1.01      0.18    -1.34    -0.67 1.00     2916     2735
## round_number    -0.04      0.01    -0.07    -0.02 1.00     2271     2736
## Condition        0.62      0.30     0.02     1.21 1.00     2793     3065
## 
## Samples were drawn using sampling(NUTS). For each parameter, Eff.Sample 
## is a crude measure of effective sample size, and Rhat is the potential 
## scale reduction factor on split chains (at convergence, Rhat = 1).

Plotting the parameters.

stanplot(b2.3)

Do those Rhat values look okay? Larger than 1.1 is deemed an issue.

rhat(b2.3) %>%
  mcmc_rhat()

How about nEff ratios? Above 0.1 is desired.

neff_ratio(b2.3) %>%
  mcmc_neff()

Only a few parameters are below the threshold.

Finally, let’s see the trace plots.

plot(b2.3, ask = F)

Looks like Stan sampled efficiently.

In this model, the fixed effect of condition is 0.62, with 95% credible intervals above 0. This implies that the participants are more likely to request from their partner in the Hidden condition. Converting to the probability scale:

post <- posterior_samples(b2.3)

visible_prob <- inv_logit_scaled(post$b_intercept)

visible_prob %>%
  median() %>%
  round(2)
## [1] 0.27
hidden_prob  <- inv_logit_scaled(post$b_intercept + post$b_Condition)

hidden_prob %>%
  median() %>%
  round(2)
## [1] 0.41
difference <- hidden_prob - visible_prob

quantile(difference,c(0.025,0.5,0.975)) %>% round(2)
##  2.5%   50% 97.5% 
##  0.00  0.14  0.28

The absolute probability difference between the conditions is +0.14 (median), with 95% CIs above 0. Participants were more likely to request from their partner in the Hidden condition.

We compute a Bayes factor for this difference between probabilities. This is the inverse of the alternative hypothesis that the two conditions are equal.

(1 / hypothesis(b2.3, "inv_logit_scaled(intercept) - inv_logit_scaled(intercept + Condition) = 0", seed = 2113)$hypothesis$Evid.Ratio) %>% round(2)
## [1] 2.38

This Bayes factor implies weak support for the hypothesis that the probabilities differ across conditions.

How much variance in the outcome does this model explain? We use the bayes_R2() function to get a Bayesian estimate of R-squared.

bayes_R2(b2.3) %>% round(2)
##    Estimate Est.Error Q2.5 Q97.5
## R2     0.27      0.01 0.24  0.29

This model explains around 27% of the variance.

4. Binomial - Probability of requesting when above the minimum threshold

Get a long-format data frame with binary request decisions over all rounds, for both conditions. If request == NA, player has died and been removed from the game, so we drop those rows. However, we also filter out rows in which the player was below the minimum survival threshold (64 cattle).

(
  d.trim <-
    d %>%
    select(ID, Group, Counterbalancing, ends_with('.player.herd_size_after_shock'), 
           ends_with('.player.request')) %>%
    gather(key, value, -ID, -Group, -Counterbalancing) %>%
    extract(key, c("Condition", "round_number", "variable"), 
            "SurvivalGame_(Hidden|Visible)[.](.|..)[.]player[.](request|herd_size_after_shock)") %>%
    spread(variable, value) %>%
    mutate(Condition    = ifelse(Condition == 'Hidden', 1, 0),
           round_number = as.integer(round_number)) %>%
    arrange(round_number, ID, Group) %>%
    drop_na() %>%
    filter(herd_size_after_shock >= 64)
)

This leaves us with 2368 request decisions.

We now fit a varying intercept and slope model, with participants nested within groups. Again, we allow the slopes for both round number and condition to vary.

b2.4 <-
  brm(data = d.trim, family = bernoulli,
      request ~ 0 + intercept + round_number + Condition
      + (0 + intercept + round_number + Condition | Group/ID),
      prior = c(prior(normal(0, 1), class = b),
                prior(student_t(3, 0, 10), class = sd),
                prior(lkj(1), class = cor)),
      iter = 2000, warmup = 1000, chains = 4, cores = 4,
      sample_prior = TRUE,
      control = list(adapt_delta = 0.99),
      seed = 2113)

b2.4 <- add_criterion(b2.4, "loo")

save(b2.4, file = "models/exp2_brmsfit4.rda")

The priors.

prior_summary(b2.4)

The results.

print(b2.4)
##  Family: bernoulli 
##   Links: mu = logit 
## Formula: request ~ 0 + intercept + round_number + Condition + (0 + intercept + round_number + Condition | Group/ID) 
##    Data: d.trim (Number of observations: 2368) 
## Samples: 4 chains, each with iter = 2000; warmup = 1000; thin = 1;
##          total post-warmup samples = 4000
## 
## Group-Level Effects: 
## ~Group (Number of levels: 40) 
##                             Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## sd(intercept)                   1.12      0.49     0.12     2.01 1.00      517      640
## sd(round_number)                0.08      0.03     0.01     0.15 1.01      397      491
## sd(Condition)                   1.93      0.76     0.31     3.34 1.01      637      677
## cor(intercept,round_number)    -0.15      0.42    -0.79     0.77 1.01      597     1505
## cor(intercept,Condition)       -0.08      0.41    -0.80     0.79 1.01      520      922
## cor(round_number,Condition)     0.33      0.38    -0.60     0.91 1.00      695      952
## 
## ~Group:ID (Number of levels: 80) 
##                             Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## sd(intercept)                   1.53      0.36     0.85     2.25 1.00      848     1412
## sd(round_number)                0.09      0.03     0.04     0.14 1.00      600      925
## sd(Condition)                   2.30      0.66     1.25     3.79 1.00      759     1842
## cor(intercept,round_number)    -0.61      0.23    -0.92    -0.02 1.00     1483     2016
## cor(intercept,Condition)       -0.25      0.25    -0.68     0.29 1.00      867     1497
## cor(round_number,Condition)     0.50      0.28    -0.14     0.93 1.01      506     1075
## 
## Population-Level Effects: 
##              Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## intercept       -1.96      0.33    -2.63    -1.33 1.00     3063     2850
## round_number    -0.10      0.03    -0.16    -0.05 1.00     2159     2447
## Condition        0.18      0.49    -0.82     1.10 1.00     3298     2919
## 
## Samples were drawn using sampling(NUTS). For each parameter, Eff.Sample 
## is a crude measure of effective sample size, and Rhat is the potential 
## scale reduction factor on split chains (at convergence, Rhat = 1).

Creating a forest plot of parameters.

stanplot(b2.4)

Viusalise Rhat values.

rhat(b2.4) %>%
  mcmc_rhat()

Viusalise nEff ratios.

neff_ratio(b2.4) %>%
  mcmc_neff()

Trace plots.

plot(b2.4, ask = F)

Looks like Stan sampled efficiently.

In this model, the fixed effect of condition is 0.18, with 95% credible intervals crossing 0. This implies that, at least when above the minimum survival threshold, the average probability of requesting did not differ between conditions.

Converting the fixed effects to the probability scale:

post <- posterior_samples(b2.4)

visible_prob <- inv_logit_scaled(post$b_intercept)

visible_prob %>%
  median() %>%
  round(2)
## [1] 0.12
hidden_prob  <- inv_logit_scaled(post$b_intercept + post$b_Condition)

hidden_prob %>%
  median() %>%
  round(2)
## [1] 0.15
difference <- hidden_prob - visible_prob

quantile(difference,c(0.025,0.5,0.975)) %>% round(2)
##  2.5%   50% 97.5% 
## -0.07  0.02  0.16

The absolute probability difference between the conditions is +0.02 (median), with 95% CIs above 0. Participants were no more likely to request from their partner in the Hidden condition.

We compute a Bayes factor for this difference between probabilities.

(1 / hypothesis(b2.4, "inv_logit_scaled(intercept) - inv_logit_scaled(intercept + Condition) = 0", seed = 2113)$hypothesis$Evid.Ratio) %>% round(2)
## [1] 0.32

This Bayes factor implies moderate support for the hypothesis that the probabilities are equal in each condition.

How much variance in the outcome does this model explain?

bayes_R2(b2.4) %>% round(2)
##    Estimate Est.Error Q2.5 Q97.5
## R2     0.43      0.02 0.39  0.46

This model explains around 43% of the variance.

4b. Binomial - Probability of requesting when above the minimum threshold (interacting with counterbalancing)

We found no effect of condition in the previous model. Is this because of order effects?

b2.4b <-
  brm(data = d.trim, family = bernoulli,
      request ~ 0 + intercept + round_number + Condition*Counterbalancing
      + (0 + intercept + round_number + Condition | Group/ID),
      prior = c(prior(normal(0, 1), class = b),
                prior(student_t(3, 0, 10), class = sd),
                prior(lkj(1), class = cor)),
      iter = 2000, warmup = 1000, chains = 4, cores = 4,
      sample_prior = TRUE,
      control = list(adapt_delta = 0.99),
      seed = 2113)

save(b2.4b, file = "models/exp2_brmsfit4b.rda")

The priors.

prior_summary(b2.4b)

The results.

print(b2.4b)
##  Family: bernoulli 
##   Links: mu = logit 
## Formula: request ~ 0 + intercept + round_number + Condition * Counterbalancing + (0 + intercept + round_number + Condition | Group/ID) 
##    Data: d.trim (Number of observations: 2368) 
## Samples: 4 chains, each with iter = 2000; warmup = 1000; thin = 1;
##          total post-warmup samples = 4000
## 
## Group-Level Effects: 
## ~Group (Number of levels: 40) 
##                             Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## sd(intercept)                   1.14      0.49     0.15     2.05 1.00      463      603
## sd(round_number)                0.08      0.03     0.01     0.14 1.01      380      583
## sd(Condition)                   1.34      0.71     0.08     2.77 1.01      345     1108
## cor(intercept,round_number)    -0.14      0.42    -0.78     0.78 1.00      591     1056
## cor(intercept,Condition)       -0.19      0.44    -0.88     0.71 1.01      872     1968
## cor(round_number,Condition)     0.31      0.42    -0.68     0.92 1.01      788     1259
## 
## ~Group:ID (Number of levels: 80) 
##                             Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## sd(intercept)                   1.53      0.38     0.83     2.28 1.00      634     1346
## sd(round_number)                0.08      0.03     0.03     0.14 1.01      467      569
## sd(Condition)                   2.14      0.55     1.18     3.29 1.01      693     1295
## cor(intercept,round_number)    -0.60      0.25    -0.92     0.02 1.00      918     1014
## cor(intercept,Condition)       -0.26      0.26    -0.69     0.33 1.01      486      935
## cor(round_number,Condition)     0.52      0.28    -0.11     0.93 1.01      399      859
## 
## Population-Level Effects: 
##                            Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## intercept                     -2.09      0.40    -2.89    -1.31 1.00     1801     2854
## round_number                  -0.10      0.03    -0.15    -0.05 1.00     2276     2460
## Condition                     -0.38      0.50    -1.41     0.56 1.00     2581     2981
## Counterbalancing               0.19      0.48    -0.79     1.10 1.00     1861     2203
## Condition:Counterbalancing     1.56      0.62     0.30     2.73 1.00     2777     2751
## 
## Samples were drawn using sampling(NUTS). For each parameter, Eff.Sample 
## is a crude measure of effective sample size, and Rhat is the potential 
## scale reduction factor on split chains (at convergence, Rhat = 1).

Creating a forest plot of parameters.

stanplot(b2.4b)

Trace plots.

plot(b2.4b, ask = F)

Looks like Stan sampled efficiently.

In this model, the interaction parameter is 1.56, with 95% CIs above zero. This implies that there is an interaction effect. Condition only has an effect on cheating when participants played the Hidden game first.

Converting the fixed effects to the probability scale:

post <- posterior_samples(b2.4b)

visible_VisibleFirst_prob <- inv_logit_scaled(post$b_intercept)

visible_VisibleFirst_prob %>%
  median() %>%
  round(2)
## [1] 0.11
hidden_VisibleFirst_prob  <- inv_logit_scaled(post$b_intercept + post$b_Condition)

hidden_VisibleFirst_prob %>%
  median() %>%
  round(2)
## [1] 0.08
difference_VisibleFirst <- hidden_VisibleFirst_prob - visible_VisibleFirst_prob

quantile(difference_VisibleFirst,c(0.025,0.5,0.975)) %>% round(2)
##  2.5%   50% 97.5% 
## -0.11 -0.03  0.06

When participants play the Visible game first, the absolute probability difference between the conditions is -0.03 (median), with 95% CIs crossing 0. Participants were no more likely to request from their partner in the Hidden condition.

We compute a Bayes factor for this difference between probabilities.

(1 / hypothesis(b2.4b, "inv_logit_scaled(intercept) - inv_logit_scaled(intercept + Condition) = 0", seed = 2113)$hypothesis$Evid.Ratio) %>% round(2)
## [1] 0.31

This Bayes factor implies moderate support for the hypothesis that the probabilities are equal in each condition.

What about for when participants play the Hidden game first?

visible_HiddenFirst_prob <- inv_logit_scaled(post$b_intercept + post$b_Counterbalancing)

visible_HiddenFirst_prob %>%
  median() %>%
  round(2)
## [1] 0.13
hidden_HiddenFirst_prob  <- inv_logit_scaled(post$b_intercept + post$b_Counterbalancing + 
                                          post$b_Condition + post$`b_Condition:Counterbalancing`)

hidden_HiddenFirst_prob %>%
  median() %>%
  round(2)
## [1] 0.34
difference_HiddenFirst <- hidden_HiddenFirst_prob - visible_HiddenFirst_prob

quantile(difference_HiddenFirst,c(0.025,0.5,0.975)) %>% round(2)
##  2.5%   50% 97.5% 
##  0.01  0.20  0.43

When participants play the Hidden game first, the absolute probability difference between the conditions is +0.20 (median), with 95% CIs above 0. Participants are more likely to request from their partner when beneath the threshold in the Hidden condition.

We compute a Bayes factor for this difference between probabilities.

(1 / hypothesis(b2.4b, "inv_logit_scaled(intercept + Counterbalancing) - inv_logit_scaled(intercept + Counterbalancing + Condition + Condition:Counterbalancing) = 0", seed = 2113)$hypothesis$Evid.Ratio) %>% round(2)
## [1] 3.05

This Bayes factor implies moderate support for the hypothesis that the probabilities differ across conditions.

How much variance in the outcome does this model explain?

bayes_R2(b2.4b) %>% round(2)
##    Estimate Est.Error Q2.5 Q97.5
## R2     0.42      0.02 0.39  0.45

This model explains around 42% of the variance.

5. Binomial - Probability of not responding to a request

Get long-format data frame with ‘received’ variable (i.e. how much a player received on any given round). We swap this around so it reflects how much the player gave to their partner (i.e. how much their partner received). We drop rows with NAs, since partners did not request help in that particular round. We then code whether the player gave nothing in response to the request (1) or gave at least one cattle (0).

(
  d.trim <-
    d %>%
    select(ID, Group, ends_with('.player.received')) %>%
    # swap every other value to get what the player GAVE, rather than what they RECEIVED
    # swapping works because each player is next to their partner in the data frame
    mutate_at(vars(ends_with(".player.received")), function(x) x[1:nrow(d) + c(1,-1)]) %>%
    gather(key, received, -ID, -Group) %>%
    extract(key, c("Condition", "round_number"), 
            "SurvivalGame_(Hidden|Visible)[.](.|..)[.]player[.]received") %>%
    rename(responded    = received) %>%
    mutate(Condition    = ifelse(Condition == 'Hidden', 1, 0),
           round_number = as.integer(round_number),
           notResponded = ifelse(responded > 0, 0, 1) %>% as.integer()) %>%
    arrange(round_number, ID, Group) %>%
    drop_na()
)

This leaves us with 716 possible responses to requests.

We then fit the varying intercept and slope model, again with participants nested within groups.

b2.5 <-
  brm(data = d.trim, family = bernoulli,
      notResponded ~ 0 + intercept + round_number + Condition
      + (0 + intercept + round_number + Condition | Group/ID),
      prior = c(prior(normal(0, 1), class = b),
                prior(student_t(3, 0, 10), class = sd),
                prior(lkj(1), class = cor)),
      iter = 2500, warmup = 1000, chains = 4, cores = 4,
      sample_prior = TRUE,
      control = list(adapt_delta = 0.99),
      seed = 2113)

save(b2.5, file = "models/exp2_brmsfit5.rda")

The priors we used.

prior_summary(b2.5)

The results.

print(b2.5)
##  Family: bernoulli 
##   Links: mu = logit 
## Formula: notResponded ~ 0 + intercept + round_number + Condition + (0 + intercept + round_number + Condition | Group/ID) 
##    Data: d.trim (Number of observations: 716) 
## Samples: 4 chains, each with iter = 2500; warmup = 1000; thin = 1;
##          total post-warmup samples = 6000
## 
## Group-Level Effects: 
## ~Group (Number of levels: 40) 
##                             Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## sd(intercept)                   1.03      0.46     0.15     2.00 1.00      992     1023
## sd(round_number)                0.09      0.04     0.01     0.18 1.01      637     1085
## sd(Condition)                   1.17      0.57     0.11     2.35 1.00      592     1054
## cor(intercept,round_number)     0.10      0.42    -0.67     0.87 1.01      814     1978
## cor(intercept,Condition)       -0.39      0.41    -0.92     0.63 1.01      720     1881
## cor(round_number,Condition)     0.36      0.41    -0.63     0.93 1.01      703     1341
## 
## ~Group:ID (Number of levels: 78) 
##                             Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## sd(intercept)                   0.54      0.41     0.02     1.53 1.00      751     1393
## sd(round_number)                0.03      0.03     0.00     0.09 1.00     1445     2133
## sd(Condition)                   0.58      0.46     0.02     1.68 1.00      708     1496
## cor(intercept,round_number)    -0.09      0.50    -0.91     0.84 1.00     2882     3262
## cor(intercept,Condition)       -0.44      0.50    -0.98     0.72 1.00     1071     2548
## cor(round_number,Condition)    -0.05      0.51    -0.90     0.87 1.00     3727     4161
## 
## Population-Level Effects: 
##              Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## intercept       -1.70      0.36    -2.45    -1.02 1.00     1879     2954
## round_number    -0.06      0.04    -0.15    -0.00 1.01     1360     2272
## Condition        0.49      0.39    -0.31     1.23 1.00     1919     2638
## 
## Samples were drawn using sampling(NUTS). For each parameter, Eff.Sample 
## is a crude measure of effective sample size, and Rhat is the potential 
## scale reduction factor on split chains (at convergence, Rhat = 1).

Plotting the parameters.

stanplot(b2.5)

Let’s look at the trace plots to make sure Stan sampled efficiently.

plot(b2.5, ask = F)

Rhat values, n_eff values, and trace plots look okay.

The fixed effect of condition is 0.49, with 95% credible intervals that cross 0. This implies that participants were no more likely to not respond to requests in the Hidden condition.

As before, we sample from the posterior and convert this difference to the absolute probability scale.

post <- posterior_samples(b2.5)

visible_prob <- inv_logit_scaled(post$b_intercept)

visible_prob %>%
  median() %>%
  round(2)
## [1] 0.16
hidden_prob  <- inv_logit_scaled(post$b_intercept + post$b_Condition)

hidden_prob %>%
  median() %>%
  round(2)
## [1] 0.23
difference <- hidden_prob - visible_prob

quantile(difference,c(0.025,0.5,0.975)) %>% round(2)
##  2.5%   50% 97.5% 
## -0.05  0.08  0.19

The probability difference between the conditions is 0.07 (median), with 95% CIs crossing 0.

We compute a Bayes factor for this difference between probabilities.

(1 / hypothesis(b2.5, "inv_logit_scaled(intercept) - inv_logit_scaled(intercept + Condition) = 0", seed = 2113)$hypothesis$Evid.Ratio) %>% round(2)
## [1] 0.67

This Bayes factor implies anecdotal support for the hypothesis that the probabilities are equal across conditions.

bayes_R2(b2.5) %>% round(2)
##    Estimate Est.Error Q2.5 Q97.5
## R2     0.24      0.03 0.18   0.3

bayes_R2() tells us that this model explains around 24% of the variance in the outcome.

6. Binomial - Probability of fulfilling a request when able

Again, the data wrangling for this model is a little trickier.

  1. First, we get a long-format data frame with (a) the player’s herd size that round, (b) how much the player received that round, and (c) how much the player requested that round.
  2. Next, we flip the latter two variables to get (b’) how much the player gave, and (c’) how much the player’s partner requested. ‘Flipping’ is possible because partners sit next to the focal player in the data frame.
  3. We drop rows with NAs, as no requesting happened this round.
  4. We keep only rows in which the partner’s request could be fulfilled without dropping the player below the minimum survival threshold (i.e. the player was able to give).
  5. Finally, we code whether the player fulfilled the request by giving what was asked or more (0) or did not fulfill the request (1).
(
  d.trim <-
    d %>%
    # 1. select herd size, received, and requested
    select(ID, Group, Counterbalancing, ends_with('.player.herd_size_after_shock'),
           ends_with('.player.received'),
           ends_with('.player.request_amount')) %>%
    # 2a. flip the variables to get what the player gave and what their partner requested
    mutate_at(vars(ends_with(".player.received"), ends_with(".player.request_amount")), 
              function(x) x[1:nrow(d) + c(1,-1)]) %>%
    gather(key, value, -ID, -Group, -Counterbalancing) %>%
    extract(key, c("Condition", "round_number", "variable"), 
            "SurvivalGame_(Hidden|Visible)[.](.|..)[.]player[.](herd_size_after_shock|received|request_amount)") %>%
    spread(variable, value) %>%
    # 2b. rename the variables to match their new meaning
    rename(gave              = received,
           partner_requested = request_amount) %>%
    mutate(Condition    = ifelse(Condition == 'Hidden', 1, 0),
           round_number = as.integer(round_number)) %>%
    arrange(round_number, ID, Group) %>%
    # 3. drop NAs, as no requesting happened
    drop_na() %>%
    # 4. was the player able to give?
    filter(herd_size_after_shock - partner_requested >= 64) %>%
    # 5. code 0 = request fulfilled, 1 = request not fulfilled
    mutate(notFulfilled = ifelse(gave >= partner_requested, 0, 1) %>% as.integer())
)

This leaves us with 473 possible response decisions in which the player was able to give their partner what they asked for. Our outcome variable is whether they fulfilled that request or not.

We now fit the varying intercept and slope model, again with participants nested within groups.

b2.6 <-
  brm(data = d.trim, family = bernoulli,
      notFulfilled ~ 0 + intercept + round_number + Condition
      + (0 + intercept + round_number + Condition | Group/ID),
      prior = c(prior(normal(0, 1), class = b),
                prior(student_t(3, 0, 10), class = sd),
                prior(lkj(1), class = cor)),
      control = list(adapt_delta = 0.99),
      sample_prior = TRUE,
      iter = 2500, warmup = 1000, chains = 4, cores = 4,
      seed = 2113)

save(b2.6, file = "models/exp2_brmsfit6.rda")

Here are the priors we used for model b2.6.

prior_summary(b2.6)

Let’s see the results.

print(b2.6)
##  Family: bernoulli 
##   Links: mu = logit 
## Formula: notFulfilled ~ 0 + intercept + round_number + Condition + (0 + intercept + round_number + Condition | Group/ID) 
##    Data: d.trim (Number of observations: 473) 
## Samples: 4 chains, each with iter = 2500; warmup = 1000; thin = 1;
##          total post-warmup samples = 6000
## 
## Group-Level Effects: 
## ~Group (Number of levels: 40) 
##                             Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## sd(intercept)                   1.20      0.57     0.13     2.34 1.00     1028      920
## sd(round_number)                0.17      0.07     0.04     0.33 1.00      889     1603
## sd(Condition)                   0.78      0.59     0.03     2.19 1.00     1314     2728
## cor(intercept,round_number)     0.22      0.39    -0.57     0.88 1.00      778     1217
## cor(intercept,Condition)       -0.03      0.49    -0.87     0.87 1.00     3004     3312
## cor(round_number,Condition)     0.22      0.47    -0.78     0.92 1.00     2221     4171
## 
## ~Group:ID (Number of levels: 73) 
##                             Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## sd(intercept)                   0.65      0.48     0.02     1.79 1.01     1125     2106
## sd(round_number)                0.08      0.05     0.00     0.19 1.00     1072     1991
## sd(Condition)                   0.77      0.59     0.03     2.20 1.00     1016     1971
## cor(intercept,round_number)    -0.18      0.49    -0.91     0.82 1.00     1601     2741
## cor(intercept,Condition)       -0.27      0.50    -0.95     0.79 1.01     1913     3600
## cor(round_number,Condition)     0.08      0.50    -0.85     0.90 1.00     2175     3743
## 
## Population-Level Effects: 
##              Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## intercept       -0.54      0.38    -1.34     0.20 1.00     2724     3783
## round_number    -0.13      0.06    -0.26    -0.04 1.00     2098     2585
## Condition        0.95      0.39     0.15     1.71 1.00     3861     3596
## 
## Samples were drawn using sampling(NUTS). For each parameter, Eff.Sample 
## is a crude measure of effective sample size, and Rhat is the potential 
## scale reduction factor on split chains (at convergence, Rhat = 1).

Plotting the parameters.

stanplot(b2.6)

Let’s look at the trace plots to make sure Stan sampled efficiently.

plot(b2.6, ask = F)

Stan sampled efficiently.

The fixed effect of condition is 0.95, with 95% credible intervals above 0. This implies that participants were more likely to not fulfill requests (when able) in the Hidden condition.

On the absolute probability scale.

post <- posterior_samples(b2.6)

visible_prob <- inv_logit_scaled(post$b_intercept)

visible_prob %>%
  median() %>%
  round(2)
## [1] 0.37
hidden_prob  <- inv_logit_scaled(post$b_intercept + post$b_Condition)

hidden_prob %>%
  median() %>%
  round(2)
## [1] 0.61
difference <- hidden_prob - visible_prob

quantile(difference,c(0.025,0.5,0.975)) %>% round(2)
##  2.5%   50% 97.5% 
##  0.04  0.23  0.40

The probability difference between the conditions is 0.23 (median), with 95% CIs above zero.

We compute a Bayes factor for this difference between probabilities.

(1 / hypothesis(b2.6, "inv_logit_scaled(intercept) - inv_logit_scaled(intercept + Condition) = 0", seed = 2113)$hypothesis$Evid.Ratio) %>% round(2)
## [1] 9.75

This Bayes factor implies moderate support for the hypothesis that the probabilities differ across conditions.

bayes_R2(b2.6) %>% round(2)
##    Estimate Est.Error Q2.5 Q97.5
## R2     0.44      0.03 0.37   0.5

This model explains around 44% of the variance in the outcome.

6b. Binomial - Probability of fulfilling request when able (interacting with counterbalancing)

As with b2.4b, we test for order effects by including an interaction.

b2.6b <-
  brm(data = d.trim, family = bernoulli,
      notFulfilled ~ 0 + intercept + round_number + Condition*Counterbalancing
      + (0 + intercept + round_number + Condition | Group/ID),
      prior = c(prior(normal(0, 1), class = b),
                prior(student_t(3, 0, 10), class = sd),
                prior(lkj(1), class = cor)),
      control = list(adapt_delta = 0.99),
      sample_prior = TRUE,
      iter = 2500, warmup = 1000, chains = 4, cores = 4,
      seed = 2113)

save(b2.6b, file = "models/exp2_brmsfit6b.rda")

Here are the priors we used.

prior_summary(b2.6b)

Let’s see the results.

print(b2.6b)
##  Family: bernoulli 
##   Links: mu = logit 
## Formula: notFulfilled ~ 0 + intercept + round_number + Condition * Counterbalancing + (0 + intercept + round_number + Condition | Group/ID) 
##    Data: d.trim (Number of observations: 473) 
## Samples: 4 chains, each with iter = 2500; warmup = 1000; thin = 1;
##          total post-warmup samples = 6000
## 
## Group-Level Effects: 
## ~Group (Number of levels: 40) 
##                             Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## sd(intercept)                   1.15      0.56     0.14     2.33 1.00     1407     1783
## sd(round_number)                0.17      0.07     0.04     0.33 1.00      845     1431
## sd(Condition)                   0.84      0.62     0.04     2.28 1.00     1249     2772
## cor(intercept,round_number)     0.26      0.38    -0.50     0.91 1.00      834     1808
## cor(intercept,Condition)       -0.08      0.49    -0.89     0.85 1.00     2254     3167
## cor(round_number,Condition)     0.23      0.46    -0.74     0.92 1.00     2665     3676
## 
## ~Group:ID (Number of levels: 73) 
##                             Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## sd(intercept)                   0.68      0.49     0.03     1.83 1.00     1380     2641
## sd(round_number)                0.08      0.05     0.00     0.20 1.00      965     2131
## sd(Condition)                   0.84      0.61     0.03     2.24 1.00     1125     2050
## cor(intercept,round_number)    -0.19      0.50    -0.92     0.84 1.00     1449     2622
## cor(intercept,Condition)       -0.30      0.49    -0.95     0.77 1.00     1877     3222
## cor(round_number,Condition)     0.08      0.50    -0.86     0.90 1.00     2420     3671
## 
## Population-Level Effects: 
##                            Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## intercept                     -0.77      0.45    -1.67     0.10 1.00     3641     3946
## round_number                  -0.13      0.05    -0.26    -0.04 1.00     2117     2799
## Condition                      0.64      0.52    -0.40     1.63 1.00     4099     3584
## Counterbalancing               0.44      0.54    -0.65     1.49 1.00     4149     4248
## Condition:Counterbalancing     0.48      0.60    -0.68     1.67 1.00     4403     4408
## 
## Samples were drawn using sampling(NUTS). For each parameter, Eff.Sample 
## is a crude measure of effective sample size, and Rhat is the potential 
## scale reduction factor on split chains (at convergence, Rhat = 1).

Plotting the parameters.

stanplot(b2.6b)

The interaction parameter has 95% CIs that cross zero, indicating no interaction effect.

Converting the fixed effects to the probability scale:

post <- posterior_samples(b2.6b)

visible_VisibleFirst_prob <- inv_logit_scaled(post$b_intercept)

visible_VisibleFirst_prob %>%
  median() %>%
  round(2)
## [1] 0.32
hidden_VisibleFirst_prob  <- inv_logit_scaled(post$b_intercept + post$b_Condition)

hidden_VisibleFirst_prob %>%
  median() %>%
  round(2)
## [1] 0.48
difference_VisibleFirst <- hidden_VisibleFirst_prob - visible_VisibleFirst_prob

quantile(difference_VisibleFirst,c(0.025,0.5,0.975)) %>% round(2)
##  2.5%   50% 97.5% 
## -0.08  0.15  0.37

When participants play the Visible game first, the absolute probability difference between the conditions is -0.03 (median), with 95% CIs crossing 0. Participants were no more likely to request from their partner in the Hidden condition.

We compute a Bayes factor for this difference between probabilities.

(1 / hypothesis(b2.6b, "inv_logit_scaled(intercept) - inv_logit_scaled(intercept + Condition) = 0", seed = 2113)$hypothesis$Evid.Ratio) %>% round(2)
## [1] 1.43

This Bayes factor implies moderate support for the hypothesis that the probabilities are equal in each condition.

What about for when participants play the Hidden game first?

visible_HiddenFirst_prob <- inv_logit_scaled(post$b_intercept + post$b_Counterbalancing)

visible_HiddenFirst_prob %>%
  median() %>%
  round(2)
## [1] 0.42
hidden_HiddenFirst_prob  <- inv_logit_scaled(post$b_intercept + post$b_Counterbalancing + 
                                          post$b_Condition + post$`b_Condition:Counterbalancing`)

hidden_HiddenFirst_prob %>%
  median() %>%
  round(2)
## [1] 0.69
difference_HiddenFirst <- hidden_HiddenFirst_prob - visible_HiddenFirst_prob

quantile(difference_HiddenFirst,c(0.025,0.5,0.975)) %>% round(2)
##  2.5%   50% 97.5% 
##  0.03  0.26  0.46

When participants play the Hidden game first, the absolute probability difference between the conditions is +0.20 (median), with 95% CIs above 0. Participants are more likely to request from their partner when beneath the threshold in the Hidden condition.

We compute a Bayes factor for this difference between probabilities.

(1 / hypothesis(b2.6b, "inv_logit_scaled(intercept + Counterbalancing) - inv_logit_scaled(intercept + Counterbalancing + Condition + Condition:Counterbalancing) = 0", seed = 2113)$hypothesis$Evid.Ratio) %>% round(2)
## [1] 7.07

This Bayes factor implies moderate support for the hypothesis that the probabilities differ across conditions.

bayes_R2(b2.6b) %>% round(2)
##    Estimate Est.Error Q2.5 Q97.5
## R2     0.44      0.03 0.38   0.5

This model explains around 44% of the variance in the outcome.

7. Plotting

Session Info

sessionInfo()
## R version 3.6.1 (2019-07-05)
## Platform: x86_64-w64-mingw32/x64 (64-bit)
## Running under: Windows 10 x64 (build 16299)
## 
## Matrix products: default
## 
## locale:
## [1] LC_COLLATE=English_New Zealand.1252  LC_CTYPE=English_New Zealand.1252    LC_MONETARY=English_New Zealand.1252
## [4] LC_NUMERIC=C                         LC_TIME=English_New Zealand.1252    
## 
## attached base packages:
## [1] grid      stats     graphics  grDevices utils     datasets  methods   base     
## 
## other attached packages:
##  [1] forcats_0.4.0   stringr_1.4.0   dplyr_0.8.3     purrr_0.3.2     readr_1.3.1     tidyr_1.0.0     tibble_2.1.3   
##  [8] tidyverse_1.2.1 ggplot2_3.2.1   cowplot_1.0.0   brms_2.10.0     Rcpp_1.0.2      bayesplot_1.7.0
## 
## loaded via a namespace (and not attached):
##  [1] nlme_3.1-140         matrixStats_0.55.0   xts_0.11-2           lubridate_1.7.4      threejs_0.3.1       
##  [6] httr_1.4.1           rstan_2.19.2         tools_3.6.1          backports_1.1.5      R6_2.4.0            
## [11] DT_0.9               lazyeval_0.2.2       colorspace_1.4-1     withr_2.1.2          tidyselect_0.2.5    
## [16] gridExtra_2.3        prettyunits_1.0.2    processx_3.4.1       Brobdingnag_1.2-6    compiler_3.6.1      
## [21] rvest_0.3.4          cli_1.1.0            xml2_1.2.2           shinyjs_1.0          labeling_0.3        
## [26] colourpicker_1.0     scales_1.0.0         dygraphs_1.1.1.6     ggridges_0.5.1       callr_3.3.2         
## [31] digest_0.6.21        StanHeaders_2.19.0   rmarkdown_1.15       base64enc_0.1-3      pkgconfig_2.0.3     
## [36] htmltools_0.3.6      readxl_1.3.1         htmlwidgets_1.3      rlang_0.4.0          rstudioapi_0.10     
## [41] shiny_1.3.2          generics_0.0.2       zoo_1.8-6            jsonlite_1.6         crosstalk_1.0.0     
## [46] gtools_3.8.1         inline_0.3.15        magrittr_1.5         loo_2.1.0            Matrix_1.2-17       
## [51] munsell_0.5.0        abind_1.4-5          lifecycle_0.1.0      stringi_1.4.3        yaml_2.2.0          
## [56] pkgbuild_1.0.5       plyr_1.8.4           parallel_3.6.1       promises_1.0.1       crayon_1.3.4        
## [61] miniUI_0.1.1.1       lattice_0.20-38      haven_2.1.1          hms_0.5.1            zeallot_0.1.0       
## [66] knitr_1.25           ps_1.3.0             pillar_1.4.2         igraph_1.2.4.1       markdown_1.1        
## [71] shinystan_2.5.0      reshape2_1.4.3       stats4_3.6.1         rstantools_2.0.0     glue_1.3.1          
## [76] evaluate_0.14        modelr_0.1.5         vctrs_0.2.0          httpuv_1.5.2         cellranger_1.1.0    
## [81] gtable_0.3.0         assertthat_0.2.1     xfun_0.9             mime_0.7             xtable_1.8-4        
## [86] broom_0.5.2          coda_0.19-3          later_0.8.0          rsconnect_0.8.15     shinythemes_1.1.2   
## [91] ellipsis_0.3.0       bridgesampling_0.7-2